Analytic solution to swing equations in power grids with ZIP load models

Objective This research pioneers a novel approach to obtain a closed-form analytic solution to the nonlinear second order differential swing equation that models power system dynamics. The distinctive element of this study is the integration of a generalized load model known as a ZIP load model (constant impedance Z, constant current I, and constant power P loads). Methods Building on previous work where an analytic solution for the swing equation was derived in a linear system with limited load types, this study introduces two fundamental novelties: 1) the innovative examination and modeling of the ZIP load model, successfully integrating constant current loads to augment constant impedance and constant power loads; 2) the unique derivation of voltage variables in relation to rotor angles employing the holomorphic embedding (HE) method and the Padé approximation. These innovations are incorporated into the swing equations to achieve an unprecedented analytical solution, thereby enhancing system dynamics. Simulations on a model system were performed to evaluate transient stability. Results The ZIP load model is ingeniously utilized to generate a linear model. A comparison of the developed load model and analytical solution with those obtained through time-domain simulation demonstrated the remarkable precision and efficacy of the proposed model across a range of IEEE model systems. Conclusion The study addresses the key challenges in power system dynamics, namely the diverse load characteristics and the time-consuming nature of time-domain simulation. Breaking new ground, this research proposes an analytical solution to the swing equation using a comprehensive ZIP model, without resorting to unphysical assumptions. The close-form solution not only assures computational efficiency but also preserves accuracy. This solution effectively estimates system dynamics following a disturbance, representing a significant advancement in the field.


I. Introduction
As societal demand for electric power grows, maintaining synchronization among synchronous generators becomes increasingly difficult. Failure to maintain synchronization might lead to fluctuations in the rotational motion of the generators over time. The swing equations [1], which are complex nonlinear second-order differential equations, govern the behavior of the synchronous generators. The problem has remained unanswered for almost a century, and research on rotor motions post-disturbance has primarily concentrated on either physically unacceptable assumptions [2][3][4][5][6][7][8][9] or time-consuming numerical simulations [10][11][12][13][14][15][16], which may not be suited for assessing real-time system stability.
An analytical solution for linear loads such as frequency-dependent loads and constant impedance loads was recently developed, allowing power grid stability to be assessed with minimal computational effort [17]. However, this approach does not take into account nonlinear loads, such as constant power loads (PQ loads). The BIG model has been proposed [18,19], but it implies load linearization and a minimal deviations from the static operating point, which may not be valid during disturbances. The deviation could lead to substantial errors in the simulation results [10,11,16]. Please see the surveys in [20,21] for further information on load modeling.
During the transients, voltages and flows are determined by power flow equations. To determine the voltages for a given set of PQ, PV, and slack buses, a power flow (PF) problem is solved. Because the PF problem is NP-hard (non-deterministic polynomial-time), finding a precise solution within an acceptable timeframe is challenging. To solve the PF problem, heuristic approaches, such as Newton-Raphson or Gauss-Seidel methods [1], are typically utilized. Each node usually has two known values and two equations for calculating the real and imaginary parts of voltages, as well as real and reactive power injections. However, a novel tensor computation-based method has been presented that does not require uniform knowledge and assigns 2N knowns from all nodes [22].
Holomorphic embedding (HE) is another method being investigated for an analytical solution to the PF problem [23,24]. HE entails equating the s term by term and extending the voltages, reciprocal voltages, and reactive power generation in terms of the variable s. The zeroth order solution, or "germ solution," is determined independently of the loading conditions, although the first and higher order terms are affected by them. The HE method is extended to numerical dynamic simulations, albeit with the strict requirements that the variables be approximated as a truncated sum of McLaurin series [25]. Because this assumption is only valid within an unknown convergence radius, the number of terms to include after truncation is determined heuristically, which potentially leads to inaccurate results if the number of terms is incorrectly assigned.
In this study, the HE method will be employed to derive voltage interdependence for various load characteristics. The following is how the paper is structured: Section II provides a load modeling approach for integrating load characteristics into swing equations in an efficient and precise manner. Section III shows how nodal voltages are linearly proportional to rotor angles (i.e., voltages at Ibus). Section IV provides the analytical solution to a differential equation established inside a validity region in order to approximate swing equations using

II.1. Constant impedance loads
Constant impedance (Z) loads have steady impedance characteristics, which means their electrical resistance remains constant regardless of the applied voltage. As a result, the current delivered to the load is directly proportional to the voltage, however the power delivered to the load increases with the square of the magnitude of the voltage. Constant impedance loads accurately depict resistance-based loads such as electric heating devices and incandescent light bulbs.
One approach to model constant impedance loads is to embed them in the Y bus matrix. The j th element in the diagonal Y add bus matrix in this method is Y add bus j; j ð Þ ¼ S 0 * j jv j j 2 [26], where S 0 * j represents the complex conjugate of the constant impedance loads at the j th node when the voltage magnitude at the node |v j | is 1.0 per unit (p.u.). When the voltage magnitude is |v j |, the actual realized power injection at Node j is S * j ¼ v j i * j � � * ¼ Y add bus j; j ð Þjv j j 2 ¼ S 0 * j jv j j 2 [26]. Consequently, the modified nodal admittance matrix is Y bus ¼ Y 0 bus þ diag S 0 * j � � , where Y 0 bus is the nodal admittance matrix without Z loads. The modified matrix Y bus is then used to calculate the network current I bus with constant impedance, which can be estimated using the formula I bus = Y bus v.

II.2. Constant current loads
Constant current loads are electrical loads that provide power proportional to the magnitude of the voltage at the node, S j ¼ S 0 j jv j j [26]. It is, however, difficult to establish the relationship between current and voltage based on the linear relationship between power and voltage magnitude. As a result, for modeling purposes, we approximated constant current loads in terms of other load types. During transients, the voltage magnitudes |v j | at the KM nodes do not deviate significantly from known values and may be expressed as jv j j ¼ jv 0 j j þ D j , where jD j j � jv 0 j j. Taylor's series expansion approximates jv j j 2 ; jv j j 2 ffi jv 0 j j 2 þ 2jv 0 j jD j , allowing constant current loads to be decomposed into constant impedance and constant power loads, The deviation of the approximation S app j =S true j is calculated as an error measurement; S app j =S true Fig 1 depicts the true delivered normalized power to constant current loads and its approximation, which contains two quadratic components in relation to voltage magnitude and do not vary with magnitude. Fig 1 shows that the approximation closely matches the true values within a range of 2-3% for a wide range of voltage magnitudes, with the approximation deviating by around 20% from the voltage magnitude at which it was made. If the voltage magnitudes exceed a specific range, the approximation may no longer be valid, and an updated approximation with the new voltage magnitude v 0 j may be necessary. However, there was little need for this in our stability studies.

II.3. Constant power loads
A constant power load is an electrical load that requires a consistent power supply to operate regardless of voltage or current fluctuations. Although the load's impedance varies, the power consumed remains constant. Electric motors, fluorescent lighting, and induction heating are examples of such loads. Voltage-dependent loads, on the other hand, are easier to support on the transmission network since they respond to variations in system voltage. Constant impedance loads, for example, fluctuate proportionally to the square of the voltage, resulting in a greater reduction when the voltage falls below 1.0 p.u. This constant power load type is observed in numerous power flow analyses, such as MATPOWER [27] and Holomorphic Embedding [23].

II.4. Concluding remark on the ZIP load models
Throughout this study, the ZIP load models are used in a variety of ways. Constant current I loads are first transformed into constant impedance and constant power loads. The extended Z and P loads are created by combining these transformed components with the existing Z and P loads. The extended constant impedance loads Z are incorporated into the nodal admittance matrix Y bus , modifying the diagonal components of the nodal admittance matrix Y 0 bus , which is generated using the branch data and shunt elements. P continues to be a load for the power flow problem or quasi-power flow problem (QPF as defined in Section III), which is discussed in the following section. As a result, the PF problem, or QPF, establishes the voltages required to fulfill the extended P loads throughout the system, using the modified Y bus .

III. Linearized interdependence between voltages
Real and reactive power injections, as well as real and imaginary voltages, are required to define the state of a node. Traditional PF problems use two known variables at each node to determine the other two, which vary depending on the bus type: 1) slack bus with known voltages; 2) PV bus with specified real power injection and voltage magnitude; and 3) PQ bus with both real and reactive power injections specified. The power balance equations at Bus i, , and the computation of the other two variables using the two known. During transients, the voltage magnitudes at Ibus stay constant [28]. The goal of QPF is to determine the voltages at KM in response to the rotor angle (equivalent to the voltage angle at Ibus). This research addresses the QPF problem, which is a form of PF problem in which the remaining two variables are calculated using power balance equations and two known variables. The QPF formulation differs slightly from standard PF problems in that it has two types of buses: reference buses with known real and imaginary voltages (all Ibus) and PQ buses.
Load characteristics are critical and this study takes ZIP loads into account. After the decomposition of I loads, QPF determines voltages throughout the system that satisfy the extended P loads with the modified nodal admittance matrix in terms of the extended Z loads. The Kronecker product is capable of efficiently solving the QPF problem [22]. Kirchhoff's cur- [1]. In this case, I 0 j is a designated current for all j in KM, serving as a placeholder to improve load model flexibility, such as shunt currents, creating a linear relationship between the voltages. Previous investigations have shown that when only Z loads are present, the linear relationship is accurate [17]. Many attempts have been made to derive a linear relationship between nodal voltages. The linear relationship has the advantage of allowing the power flow equations to be concisely expressed in a smaller number of nodes, which is why it is sometimes referred as a network reduction. Section III.1. illustrates conventional techniques, and Section III.2. presents the proposed approach based on Holomorphic Embedding.

III.1. Conventional methods to represent v KM in terms of v I
The Ward reduction method [29] is employed when the voltage magnitudes at KM are almost constant,. The current at KM can be expressed as a linear combination of the voltages at KM and Ibus, This method can potentially reduce the complexity of the power flow problem while increasing computational efficiency. Voltage sensitivities close to the original dispatch can also be investigated: (v 0 ) = 0, must be maintained while voltages fluctuate. As a result, the first-term Taylor series expansion provides a linear relationship between Ibus and other buses' voltages: Because both methods investigate physical power flow models, a reduced network with fewer nodes can be constructed. However, if the voltages at Ibus (i.e., rotor angle) differ from the voltages used for linearization, the errors in the voltage estimates can be substantial.

III.2. Proposed method based on the holomorphic embedding
The holomorphic embedding method is an innovative approach to power flow analysis that tries to discover voltage solutions analytically rather than iteratively [23]. In HE's load modeling, constant impedance (Z) load models are represented as I = Y bus v + I 0 where I represents current, Y bus is the nodal admittance matrix, v represent voltages, and I 0 is the designated current. The Y bus matrix has been altered and divided into two parts: Y tr and Y sh . Y tr has zero row sums meaning that its row sums are all zeros. In contrast, Y sh is a diagonal matrix made up of nodal shunt elements. If Y sh is zero, the germ solution serves no flows or loads in the system. Voltage magnitudes and angles in the germ solution are homogeneous, with a typical magnitude of unity and an angle of zero. It is possible, however, to select a reference voltage magnitude other than unity, which is denoted by E ref .
Depending on the bus types specified in [17], the power balance equation is (1A) for Bus i 2 Ibus and (1B) for Bus j 2 KM: Here, E i os the voltage magnitude at Bus i; p i and q i are the real and reactive power injections at Bus i, while p j and q j are the real and reactive power injections at Bus j; and w j is the reciprocal voltages at Bus j, i.e., w j v j = 1, respectively.
The formulae for Y tr v and Y sh v, given the partitioned into submatrices Y tr and Y sh are as fol- . ., n-1 where A is a row vector and a is a scalar. We can write w * j n ½ � as: Where M w * j n (2A) is quadratic in Ibus voltages, but it can be linearized around known bus voltage values as follows: Where and Where e T j B v KM nþ1 Where ½ � are n KM � n I and n KM × 1, respectively. Summation over n gives: x KM y KM Because adding an infinite number of terms is infeasible, we seek the Padé approximation [24], which was inspired by [23]: In this study, we set l equal to m. The first linear system is constructed using the coefficients of s k , k = 0, . . ., l: For For s l : The second linear system is given by the coefficients of s k , k = l+1, . . ., l+m: For s l+2 : For s l+m : It is worth noticing that the linear matrix equations in (3E)-(3G) are derived in such a way that the superposition s becomes zero at any value for x I and y I . To find a linear equation, use (3E)-(3G): As a result, the vector formed of the coefficients in the denominator series of the Padé approximation resides in the null space of the W b matrix. Because the denominator function of the Padé approximation (3A) cannot be zero, the vector must be a nontrivial solution. To explore the matrix's null space, the singular value decomposition of The coefficients in the nominator series are found utilizing the values for b v j k ½ � from (3B)-(3D), and (4B). The Padé approximation yields a linear relationship between v I and v KM at s = 1: Where The power balance equation for Ibus is linearizing the right-hand side of the power balance equation at Ibus using (5) and As with (3A)-(3I), the Padé approximation is applied to produce the linear expression for real and reactive power generation at Ibus: The definitions of ξ i (� q i x i ) and ξ i (� q i y i ) reveal the following: Where H I Q k , the Padé approximation offers the linear formula for ξ I and ξ I at Ibus: The linearization results are summarized below:

III.3. Error analysis
The number of terms in the nominator and denominator series, denoted by l and m, respectively, influences the accuracy of the Padé approximation. While not as robust as the McLaurin series, the amount of terms utilized has an effect on the precision of the Padé approximation [30]. The holomorphic embedding procedure limits the series to a limited number of terms, resulting in errors in the matrices and vectors due to inaccuracies in the data and observations. The errors are caused by two factors: deviation from the voltage values at Ibus during the linearization process and truncation error due to the finite nature of l. The matrix and currents obtained using the Padé approximation are represented by C l and c l , respectively. In contrast, C 1 and ψ 1 denote the matrices and currents acquired by employing an infinite number of terms. Denote w true v true as the solutions with and without errors in (9): The difference between solutions is bounded by [31]: Where � ¼ max ; and (10) is an equality constrained least squares problem with a fixed voltage reference angle, as given in [31]:

IV. Analytical solution to swing equations
When there is a disturbance in a power system, the swing equation governs the system dynam- where M j and D j are the inertia and damping coefficients of the j th generator, respectively; δ j , p mech j , and p elec j are the rotor angle, mechanical power input, and electrical power output of the j th generator, respectively. Kirchhoff's laws govern the output of electrical power. The coordinate transformation of the swing equation to [17] yields: Using (9), (11) becomes: Where It should be noted that the factors in square brackets in (12) and (13) are related to the admittance matrix and are normalized to the inertia of the corresponding generators. (12) and (13) are both simplified [17]:  (14), like the one in [17], is constrained by two parameters, O1 and O2. Because of these constraints, (14) is equivalent to the swing Eq (10). O1 denotes all of the voltages at Ibus, whereas O2 denotes the boundary condition where the linear approximation (9) is still valid. (16) depicts the combination of O1 and O2:

O1
: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 À 1 E j ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi The validity region is defined by (16), as indicated in [17]. When (16) is satisfied, (14) remains within a distance of ε from the solution to the swing equation. As a result, (14) is also ε-close to the true solution of the swing equation, and the error is bounded by k (17) ε, where k (17) is a constant as proved in [17]. According to [32], velocity control in the form of x j dx j dt þ y j dy j dt ¼ 0 is more effective in maintaining O1 than location control. The following is the solution to differential Eq (14): Where, as demonstrated in [17], Ψ k is the matrix form of the k th vector spanning the null space The rotor angle δ I is determined as the imaginary component of log(x I + jy I ). Because the reference angle bus is fixed while calculating the sensitivities in (9), the rotor angles are relative to the reference angle node ref. To put it another way, the estimated rotor angles are Using the definition of δ COI , it is possible to determine: (18) yields (19), which alters the reference angle: δ ij is the Kronecker delta [33]. The rotor angles relative to the moment of inertia can be determined directly using the solution in (17) since the linear coefficients, d ij À M i M tot , are independent of the reference bus. As a result, the linear approximation using HE in (9) remains valid, despite the fact that it is predicted on a specific choice of the reference angle bus.

V. Beyond the validity region
Denote f and g as (14) and (15), respectively: F(w I , ω I , ρ I ) is formed as follows: According to the definition of the validity region in [17], the value of F is non-zero at the last point on the boundary. Iteratively, the Newton-Raphson method is used to obtain a consistent initial point for w I ; o I ¼ dw I dt À � , and r I ¼ do I dt À � . The starting point needs to be compatible with the differential equation (14). Kirchhoff's laws are utilized to compute the initial power injection.
The solution in (17) is valid if the power injection is close enough to the estimated value from (9) and (17) as well as the update for l. Otherwise, the linear sensitivities in (9) must be recalculated to determine matrices H and L, and the solution in (17) must be recalculated, resulting in greater complexity.

VI.1. Exact approach: Time-Domain Simulation (TDS)
The , are then calculated. . This entails calculating the power generated at each node, which is then used to determine the second-order derivatives of rotor angle, dω i , respectively. Solving the QPF problem determines power generation.

Return to
Step 1 for the next iteration based on the revised voltage angles at Ibus δ i 3) .
In terms of node categorization, the QPF problem at Step 1 differs from conventional PF problems. There are reference nodes, PV nodes, and PQ nodes in conventional PF problems. In the QPF problem, the voltage magnitudes and angles are known at all Ibus while the real and reactive power injections are known at the remaining nodes, which are referred to as KM. Normally, the voltage angle at the reference angle bus is set to zero.
The first step of the QPF problem has 2N equations and 2N unknowns that can be utilized to compute the vector [v x ; v y ; p inj ; q inj ] 4N , where v x and v y are the real and the imaginary components of the voltage vector, and p inj and q inj are the real and the reactive power injections, respectively. Despite the fact that Step 1 differs from the conventional PF problem in terms of node classification, it can still be solved using a method such as the Kronecker product [22]. The Kronecker product is a matrix operation that can convert power flow equations into matrices, making numerical computation easier to solve the equations. The 2N equations and 2N unknowns in Step 1's QPF problem can be written in matrix form using the Kronecker product and solved using a variety of numerical approaches. However, the computational costs of numerically solving (23) increase considerably with each iteration. For example, 10 4 instances of (23) must be solved to complete a 10-second event with a time step of 10 −3 seconds.

VI.2. Heuristic Truncated McLaurin Series (HTMS) approach
Section V.1 emphasizes that the QPF problem is structured in the same way as conventional PF problems. To investigate this similarity, the HE [25] is employed for numerical simulations. One may begin by assuming that the time-dependent rotor angles for both δ j and ω j as a McLaurin series near t = 0: The series is then truncated to the first m terms.
With (24), one can deduce expressions for do j =dt ¼ P mÀ 1 k¼0 k þ 1 ð Þo j k þ 1 ½ �t k and Comparing the τ th power of t in (24) and the definitions of δ and ω leads to (25) that establishes the relationship between two consecutive terms, τ and τ + 1: At each term, the rotor angles are computed using (25), and two nodal variables, either nodal voltages or power injections, are determined based on the bus type, resulting in a QPF problem.
The numerical simulation process begins with a simulation of the system at τ = 0, ω j [0] and δ j [0] are calculated. The HE algorithm uses these values to find an analytic solution to identify voltages and power injections p elec j 0 ½ � at τ = 0. The values for p elec j 0 ½ � and (25) allow us to calculate ω j [1] and δ j [1] for the next HE at τ = 1, and so on.
In (24), the state and algebraic variables u(t) are assumed in a truncated McLaurin series form; u t ð Þ ¼ P m k¼0 u k ½ �t k . The series' convergence radius is: Because the series used in HE are numerically evaluated, there is no theoretical limit to the value of k min necessary to satisfy (27A) as long as k is greater than k min . As a result, the number of truncations is heuristically chosen, despite the fact that the number is important to the accuracy of the approximation. For example, the state and algebraic variables x and y are written as P 4NI j¼1 c j exp l Re j þ jl Re j � � t h i [17]. The resulting rotor angle (voltage angles at Ibus) is given by . The time dependence is clearly not polynomial.
Considering two extreme cases, rapidly diverging and rapidly converging, we observe that the function is more accurately approximated as the number of terms m in the McLaurin series increases. See Fig 2 that depicts the truncated McLaurin series of exp(t)cos(3t) (exponentially diverging case) and exp(-t)cos(3t) (exponentially converging case). The convergence region for the truncated McLaurin series for δ j , ω j , x i , and y i should be the minimum of all k min values defined in (26) at the same time, as given by (27) to ensure proper convergence for all variables:

for all i; j; and m ð27Þ
Because the HE method is based on numerical term-by-term computation, deriving the analytical forms of To the best of our knowledge, predicting t (corresponding to m) in advance to satisfy (27) is thus impossible. When the truncation error is too significant, it can be decreased by either 1) increasing the number of terms (i.e., increasing m) or 2) using the Taylor series expansion when the truncation error becomes too significant, such as To keep the error within an acceptable range, m may need to increase exponentially as the simulation's termination time t te increases. As a result, the second approach is employed in this study. Because this method does not involve updating the state and algebraic variables at each iteration, its computation cost is lower than that of traditional time-domain simulation methods. The accuracy of this approach, however, may be jeopardized if the system evolves to a new equilibrium in which the McLaurin series no longer accurately approximates the system dynamics. In such instances as shown in [25], this approach may produce imprecise numerical results or incur a high computational cost.

VII.1. Load modeling ZIP loads to ZP loads
Constant current loads are split into the constant impedance and the constant power loads using the formula S app . Power flow studies are performed on all the systems in the MATPOWER [27] to identify voltage caused errors caused by the approximation of the constant current loads. The errors remain less than 1% for all the systems. In the simulations, loads are equally distributed as Z, I, and P loads, i.e., one third of loads are Z loads, I loads, and P loads each at all the PQ nodes. Fig 3 depicts the typical examples from the selected systems.
If the voltage magnitudes do not deviate more than 10%, relative errors in the loads, as shown in Fig 1, remain below 1%. As a result, the resulting voltage errors are also less than 1%. The numerical error in approximating loads from constant current loads, as depicted in Fig 1, falls within a range similar to that of the demand forecast. The voltage difference is normally minor as long as demand forecast is close to the actual demand. In general, the resulting voltages are determined by how the load is distributed across the network. If the input is close to a certain value, the system's output will remain close to the expected output, as illustrated in

VII.2. Voltage estimate using holomorphic embedding
Section III discusses several circuit analysis methods for simplifying complex electrical networks. These methods represent an arbitrary two-port network is represented by these methods as a single voltage source in series with an impedance. By utilizing circuit analysis techniques, they enable determination of voltage and current behavior at various points in the network. The resulting equivalent circuit provides a mean to model the original circuit's behavior in terms of voltage, which proves valuable in reducing electrical circuit complexity and enhancing the understanding of their behavior. The networks considered for numerical simulation are the same as in Section VII.1, with a typical example displayed in Fig 4. All tested networks perform similarly.
The proposed method clearly provides an accurate estimate of the voltages at KM while covering a wide range of voltages at Ibus. Considering the accuracy of the estimate and the formulation of all nodal voltages and branch flows are expressed in terms of voltages at Ibus, it is possible to construct an efficient, accurate, and highly adaptable "reduced" network that maps the power flows over the original network. However, the challenge of developing a network corresponding to the highly asymmetric nodal admittance matrix appearing in I = Yv poses a potential barrier to creating a "reduced" network.

VII.3. Simulation environments
Simulations were conducted on various IEEE model systems (including the IEEE 4-, 9-, 14-, 30-, and 39-bus systems). However, for purpose of visual presentation and comparison with the results from [17], the simulation results for the IEEE 9-bus system are discussed in detail (Fig 5). The 9-bus system is expanded to 12 buses by adding 3 Ibus (Bus 10, 11, and 12) to represent the internal generators located at Bus 1, 2, and 3, respectively.
In one particular case, a line failure is assumed to occur near Bus 7 on the line connecting Bus 7 and Bus 8, which consist of two separate circuits. The fault is detected immediately after the disturbance, and the circuit breaker is opened to isolate the disturbance, indicating that the fault is cleared. Once the fault is cleared, power is restored to the remaining circuit. The stability is assessed using a coupled oscillator model [34], direct methods [2-7, 35, 36], numerical TDS, and the proposed method. Other outages, such as various load losses at Bus 8, are modelled and analyzed for their effects on system stability simulated and evaluated using the same methodologies.
When the voltage magnitudes at all nodes deviate from the pre-fault state along the on-fault trajectory of the line outage before the fault is cleared, the errors in load modeling, as well as the errors from the linearization in the holomorphic embedding, become considerable. In contrast to the errors discussed above, the inaccuracy induced by assuming that all loads are constant impedance loads is rather minor. Consequently, utilizing the constant power load type during the on-fault trajectory is feasible. However, because the fault-clearing process is quick, precise TDS does not necessitate a high computation cost. When the loads are not entirely constant impedance loads, the on-fault trajectory is identified using TDS simulation in this study.
It is worth noting that the HTMS method converges within a time radius of around 10 −2 -10 −1 seconds, with 20-30 terms yielding accurate results when compared to TDS results, with a 5% tolerance between the two components. Each HE procedure at involves a matrix inver- All numerical calculations were performed on a Mac Pro equipped with two 2.93 GHz 6-core Xeon processors. The P load values are adjusted to match the Z load values in the postfault trajectory presentedin [17], i.e., s j Zload ¼ jv j j 2 s j base and s j Pload ¼ s j base , where v j denotes the voltage at Bus j at the beginning of the post-fault trajectory in [17]. Consequently, the stability assessment results in Ref [17] for the coupled oscillator and the direct method remain unchanged. Table 1 only shows the evaluation results for TDS and the proposed method.

VII.4. Line outage between Bus 5 and Bus 7 near Bus 7
The line outage results will serve as an example. As shown in Table 1, the TDS approach is employed to calculate the on-fault trajectory, which increases the computation time of the  proposed method. In contrast to Z loads where the loads gradually adjust based on voltage magnitudes, the rotor motions are constrained to support constant power loads. Because of these constraints imposed on P loads, the real and imaginary components of the rotor motion (the time derivatives of the voltages) deviate significantly from zero. Voltages at Ibus (rotor) with P loads, as a result, cross the validity region boundaries more frequently than voltages at Ibus (rotor) with Z loads do to satisfy (16). The validity region is a key element in the proposed method since it assures that the resulting solution obtained is sufficiently close to the time series provided by TDS. In the case of the line outage simulation, the validity region boundary 1 À 1 E j ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi  Fig 7(C) compares the blue and red lines showing that when the value of ε is set to 1% the approximation employed to simplify O1 holds well. The 1-norm is used to derive a useful expression for the validity boundary based on the inequality of norms [31], k�k 1 � k�k 2 � k�k 1 .
The time where the time series intersects the validity boundary has an analytic form, and the consistent initialization yields ρ I , ω I , and v I . Using the consistent initial values, the swing equations and corresponding solutions are identified for the next validity region, ensuring that the analytical solution and numerical computation remain sufficiently close to evaluate the stability over the entire range. The stability assessment is as follows. With a known analytical solution valid within a boundary, one can assess if the system is stable or becomes unstable within the boundary. If the system becomes unstable, the assessment determines that the system is unstable. Otherwise, the assessment is continued to the next boundary. If the assessment is not completed by the time frame that a control mechanism, such as governor control, is engaged, the system is at least operationally stable [17]. Only the value for β in (17) should be modified to match the specified values for analytic solution. The computation costs associated with the update of β are minimal because the evaluation of β is a least square problem with a fixed matrix, [ψ 1 � � �ψ 4NI ].
It is important to note that the computation time in Table 1 includes all procedures until stability is determined, which entails approximately 10 times of consistent initialization problems until t = 3 seconds. This considerably increases the computation cost. The proposed approach, however, has a substantially reduced computation cost when compared to TDS. A greater value for ε in (16) can be utilized to further improve computational efficiency at the expense of compromising analytical solution accuracy. Table 1 shows the computation times for TDS and the proposed method until stability is assessed (t = 1 second). When ε is set to 10%, the initialization problem is only solved twice, while the stability assessment remains unaltered, resulting in a reduction in computation time from 3.71 seconds to 1.22 seconds.

VII.5. Voltage linearization using holomorphic embedding
The HE linearization method relies on a specified reference angle, typically zero, at an arbitrarily chosen node (in this study, Bus 1). To achieve this, we rotate the voltages until the angle at the reference node is zero. Following the fault clearance (t = 0 sec), the rotor angles undergo step changes, as shown in Fig 8(A) and 8(B). These changes result from the voltage angle shift at the reference angle node, not from alteration in the physical rotor angles. The rotor angles with respect to that of COI offsets the step change as illustrated in (19) (See Fig 8(B)). TDS will https://doi.org/10.1371/journal.pone.0286600.g007 continue until the system establishes a new equilibrium or becomes unstable. Each TDS step necessitates solving 2 QPF problems-one at the predictor stage and the other at the corrector stage, for a total of 6,000 QPF are solved during the 3-second at a frequency of 10 −3 sec. As the proposed method offers an analytic expression for the variables over time during the post-fault period, it is possible to calculate errors. For example, the leading terms in the linearization error based on HE are of the order ϑ(NKM 2 ) [31] if Y KMKM tr À � À 1 is stored. Fig 9 illustrates nodal voltages in the IEEE 9-bus system, where the imaginary components of the voltages at selected nodes (Bus 1, 4, 5, and 6) are non-positive. This observation implies that the behavior of the nodal voltages may not be uniform across the system, indicating that there may be an underlying structure that warrants further investigation. To investigate this further, the sensitivity matrix H KM , as given in (9)  is defined as θ ij . Two spaces are considered nearly linearly dependent when the angle between them is small [37,38] (See Table 2 for the values for θ ij ). Partitioning the buses based on the behavior exhibited in Fig 9 may be difficult because θ ij represents the angle between two subspaces spanning 2NI dimensions. To cluster the buses, the k-means clustering algorithm [39,40] is employed, with the θ ij values. According to k-means clustering results, Buses 1, 4, 5, and 6 constitute one group, while the remaining buses form another. Fig 5 displays the locations of these groups in the IEEE 9-bus system, with the lower half of the system forming one group. Following a line outage, the path between the two groups becomes less efficient for power transfer, leading to a degree of system partitioning. This partitioning impact is reflected in the sensitivity matrix H KM , which can provide valuable insights for analyzing the system's stability and vulnerability during disturbances. Further studies can focus on understanding and mitigating the impacts of these partitioning effects on power system performance and reliability. , and constant power loads (P), using the proposed method, which yields results that are visually indistinguishable from those obtained using the TDS method. This implies that in this study, the two major error sources of the proposed analytical solution to swing equations-I load modeling in terms of Z and P loads, and the condition in (16)-are minimal. Fig 10(A) initially displays load characteristics before 0.5 seconds, which become evident after 2 seconds (for further details, see Fig 10(B)-10(D)). According to Section II.3, higher nodal voltages at a node necessitate a larger power for P loads. This is reflected in Lw I in (14), resulting in a larger L matrix.

VII.6. The effects of load characteristics on power system dynamics
The analysis of the eigenvalues and coefficients b * k C k from (17) reveals that only two coefficients associated with the generator at Bus 3 have significant values (see second row of Table 3). Their corresponding eigenvalues are listed in the third and the fourth rows of Table 3. The first set of eigenvalues have larger imaginary components, suggesting that the corresponding functions u in (17) oscillate more quickly than those related to the second set. Therefore, the smaller coefficients of the second set of imaginary eigenvalues cause the modes' linear sum to oscillate even faster. This explains why the periods of the rotor dynamics decreases with I and P loads, potentially due to the damping effect of the constant power loads on the system.
To simplify our physical explanation, let's disregard l in (14) and suppose that all coefficients in (14) are scalar. With this adjustment, (14) then represents a spring-mass-dashpot system, where a spring either push or pull the mass illustrated in Fig 11. The equation of the motion for this system is: m d 2 dt 2 x þ c d dt x þ kx ¼ 0 [33] where m denotes the mass, c is the dashpot constant, and k represents the spring constant. The solution to this equation is: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi c 2 À 4mk p � � and l 2 ¼ 1 2m À c À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi The values for x 1 and x 2 are selected to satisfy the system's initial and boundary conditions. When c 2 is less than 4mk, the system oscillates underdamped.     7. In case where k = 1 (represented by the blue line), the system is overdamped since c 2 > 4mk, and the mass shifts to the next stable point. The red and green lines display changes in the mass location for lower and higher k values respectively. Importantly, Fig 12 shows that as k value increases, the oscillation period decreases. Fig 10 shows that the oscillation periods τ change with the types of loads, τ Zload > τ Iload > τ Pload . This suggests that in the spring-mass-dashpot system, the loads function similarly to a spring, and constant power loads exert the strongest restoring force.

VIII. Conclusions
The author of this paper proposes a novel approach for evaluating the transient stability of power grids by deriving an analytic solution to the swing equation that can accommodate all ZIP models while complying to Kirchhoff's laws. The swing equation is transformed into a differential algebraic equation, and that the ZIP load models are analytically integrated using the holomorphic embedding method. As demonstrated in the paper's example, the proposed solution enables an accurate and efficient assessment of stability.
The proposed assessment tool is fundamentally different from traditional methods such as the coupled oscillator model and the direct method because it avoids making physically inadmissible assumptions. Additionally, it distinguishes itself from time-domain simulation methods by offering relatively low computation costs. These advantages make the proposed tool an attractive option for real-time transient stability assessment in power systems, which is crucial for maintaining grid reliability and security.
The potential applications of this novel method include real-time monitoring and control of power systems, as well as grid planning and optimization. By providing a more accurate and computationally efficient means of assessing transient stability, the proposed approach can help grid operators make better-informed decisions and respond more quickly to potential threats, ultimately enhancing the overall performance and resilience of power systems.